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Abstract 

In this work, the dependence of reversal potentials and zero-current fluxes on diffusion 
coefficients are examined for ionic flows through membrane channels. The study is con¬ 
ducted for the setup of a simple structure defined by the profile of permanent charges with 
two mobile ion species, one positively charged (cation) and one negatively charged (anion). 
Numerical observations are obtained from analytical results established using geometric sin¬ 
gular perturbation analysis of classical Poisson-Nernst-Planck models. For 1:1 ionic mixtures 
with arbitrary diffusion coefficients, Mofidi and Liu \ arXiv:1909.01192 ] conducted a rigorous 
mathematical analysis and derived an equation for reversal potentials that, in its particu¬ 
lar case, can be compared to Goldman-Hodgkin-Katz equation. We summarize and extend 
these results with numerical observations for biological relevant situations. The numerical 
investigations on profiles of the electrochemical potentials, ion concentrations, and electrical 
potential across ion channels are also presented for the zero-current case. Moreover, the be¬ 
havior of current and fluxes with respect to voltages and permanent charges are investigated. 

In the opinion of the authors, many results in the paper are not intuitive, and it is difficult, 
if not impossible, to see all cases without investigations of this type. 
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1 Introduction. 

Ion channels are proteins found in cell membranes that create tiny openings in the membrane 
to allow cells to communicate with each other and with the outside to transform signals and to 
conduct tasks together Him They have an aqueous pore, that becomes accessible to ions after 
a change in the protein structure that makes ion channels open. Ion channels permit the selective 
passage of charged ions formed from dissolved salts, including sodium, potassium, calcium, and 
chloride ions that carry electrical current in and out of the cell. 

To unravel how ion channels operate, one needs to understand the physical structure of ion 
channels, which is defined by the channel shape and the permanent charge. The shape of a 
typical ion channel is often approximated as a cylindrical-like domain. Within a large class of 
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ion channels, amino acid side chains are distributed mainly over a “short” and “narrow” portion 
of the channel, with acidic side chains contributing permanent negative charges and basic side 
chains contributing permanent positive charges. The spatial distribution of side chains in a 
specific channel defines the permanent charge of the channel. 

The most basic function of ion channels is to regulate the permeability of membranes for a 
given species of ions and to select the types of ions and to facilitate and modulate the diffusion 
of ions across cell membranes. At present, these permeation and selectivity properties of ion 
channels are usually determined from the current-voltage (I-V) relations measured experimen¬ 
tally ( |13L 21]). Individual fluxes carry more information than the current, but it is expensive 
and challenging to measure them ([25| [29]) - The I-V relation defines the function of the channel 
structure, namely the ionic transport through the channel. That transport depends on driving 
forces expressed mathematically as boundary conditions. The multi-scale feature of the problem 
with multiple physical parameters allows the system to have great flexibility and to exhibit vi¬ 
brant phenomena/behaviors - a great advantage of “natural devices” ([12]). On the other hand, 
the same multi-scale feature with multiple physical parameters presents an extremely challenging 
task for anyone to extract meaningful information from experimental data, also given the fact 
that the internal dynamics cannot be measured with present techniques. The general inverse 
problem is challenging, although specific inverse problems have been successfully solved with 
surprisingly little difficulty using standard methods and software packages ([5]). 

To understand the importance of the relation of current and permanent charges, i.e., I-Q, 
we point out that the role of permanent charges in ionic channels is similar to the role of doping 
profiles in semiconductor devices. Semiconductor devices are similar to ionic channels in the way 
that they both use atomic-scale structures to control macroscopic flows from one reservoir to 
another. Ions move much as quasi-particles move in semiconductors. Semiconductor technology 
controls the migration and diffusion of quasi-particles of charge in transistors and integrated 
circuits. Doping is the process of adding impurities into intrinsic semiconductors to modulate 
its electrical, optical, and structural properties ([T7, 153]). 

Ion channels are almost always passive and do not require a source of chemical energy 
(e.g., ATP hydrolysis) in order to operate. Rather they allow ions to flow passively driven by a 
combination of the transmembrane electrical potential and the ion concentration gradient across 
the membrane. For fixed other physical quantities, the total current X = X(V, Q) depends on 
the transmembrane potential V and the permanent charge Q. For fixed Q, a reversal potential 
V = Vrev(Q) is a transmembrane potential that produces zero current Z(V re v(Q), Q) = 0. 
Similarly, for fixed transmembrane potential V, a reversal permanent charge Q = Q re v(V) is a 
permanent charge that produces zero current X(V, Qreti(V)) = 0. 

The Goldman-Hodgkin-Katz (GHK) equation for reversal potentials involving multiple ion 
species ( [22l 124] ) is used to determine the reversal potential across ion channels. The GHK 
equation is an extension of the Nernst equation - the latter is for one ion species. The classical 
derivations were based on the incorrect assumption that the electric potential <f>(A) is linear in X 
- the coordinate along the longitude of the channel. This assumption is particularly unfortunate 
because it is the change in the shape of the electrical potential that is responsible for so much 
of the fascinating behaviors transistors or ionic systems (dnum sa esi eoi Eg). There was no 
substitute for their equations until authors of da S3] recently offered equations derived from 
self-consistent PNP systems, to be the best of our knowledge. 

There have been some achievements recently in analyzing the PNP models for ionic flows 
through ion channels ([HI EH] CHI El] ES ED], etc.). Although mathematical analysis plays 
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a powerful and unique role to explain mechanisms of observed biological phenomena and to 
discover new phenomena, numerical results are needed to fit actual experimental data and 
study cases where analytical solutions do not exist. Furthermore, numerical observations give 
clues for more theoretical investigations. Indeed, numerical and analytical studies are linked; 
any progress in one catalyzes work in the other. 

The numerical results, throughout the paper, are gained from the algebraic systems (2.1), 


(2.7), (2.8) and (3.3), obtained from reduced matching systems of analytical results in 
and m- The nonlinear algebraic systems are then solved by the Matlab® function fsolve 
that uses the trust-region dogleg algorithm. The trust-region algorithm is a subspace trust- 
region method and is based on the interior-reflective Newton method described in [9]. Our 
numerical results indicate that current-volt age and current-permanent charge and even zero- 
current relations depend on a rich interplay of boundary conditions and the channel geometry 
arising from the mathematical properties analyzed in ( [14L [3T1 f43l 154] ). 

This paper is organized as follows. The classical PNP model is provided for ionic flows 
in Section 1.1 to prepare the stage for investigations on the next sections. In Section [2j we 
study zero current problems to investigate zero-current fluxes, and reversal potentials V rev . In 
particular, we compare a special case of the reversal potential with the GHK equation. Some 
other numerical observations are also provided to study profiles of relevant physical quantities 
In Section [3j we first recall the analytical results in [H] when diffusion constants 


in Section 2.4 


are also involved. Then numerical observations are provided to examine behaviors of current, 
voltage, and permanent charge with respect to each other in some general cases. Some concluding 
remarks are provided in Section [4] 


1.1 Poisson-Nernst-Planck models for ionic flows. 

The PNP system of equations has been analyzed mathematically to some extent, but the equa¬ 
tions have been simulated and computed to a much larger extent (0 0 [8, 26} [28]). One can 
see from these simulations that macroscopic reservoirs must be included in the mathematical 
formulation to describe the actual behavior of channels (120, IQ). For an ionic mixture of n ion 
species, PNP type model is, for fc = 1,2, 


Poisson: 


Nernst-Planck: 


V • (e r (X) £o V$) = -e 0 ( ^ z s C s + Q(X)), 

5—1 

dtCk + V • Jk = 0, — Jk = — —T>k{X)CkS7Hk, 

k B T 


( 1 . 1 ) 


where X 6 P with P being a three-dimensional cylindrical-like domain representing the chan¬ 
nel of length L (nm = L x 10 _9 m), Q(X) is the permanent charge density of the channel 
(with unit 1M = 1 m ° l = 10 3 ^^ ), £ r (X) is the relative dielectric coefficient (with unit 1), 
£o(~ 8.854187817620 x 10 -12 fm -1 ) is the vacuum permittivity, eo(~ 1.60217646 x 10“ 19 
C =coulomb) is the elementary charge, k B (= 1.380648813 x 10~ 23 J/\ _1 ) is the Boltzmann 
constant, T is the absolute temperature (T ~ 273.16 K =kelvin, for water); is the electric 
potential (with the unit V = Volt = JC _1 ), and, for the k -th ion species, Ck is the concentra¬ 
tion (with unit M), Zk is the valence (the number of charges per particle with unit 1), ^k is the 
electrochemical potential depending on 4> and Ck (with unit J = CV). The flux density Jk[X) 
(with unit Hz = f) is the number of particles across each cross-section in per unit time, Vk(X) 
is the diffusion coefficient (with unit m 2 /s), and n is the number of distinct types of ion species 
(with unit 1). 
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Ion channels have narrow cross-sections relative to their lengths. Therefore, three-dimensional 
PNP type models can be reduced to quasi-one-dimensional models. The authors of [H) first of¬ 
fered a reduced form, and for a particular case, the reduction is precisely verified in [37] . The 
quasi-one-dimensional steady-state PNP type is, for k = 1,2, ...,n, 


1 


d 


A(X ) dX 


e r (X)e 0 A(X)^j = - e 0 


J2z s C s + Q(X) 


\S= 1 


dJk 

~dX 


= o, -Jk = 


X Vt(x)A{x)Ct ^ 


( 1 . 2 ) 


where A(X) is the area of cross-section of the channel over location X. Equipped with system 

,n, 


(1.2), we impose the following boundary conditions, for k = 1,2, 

$(0) = V, C'fc(O) = L k > 0; 4>(L) = 0, C k (L) = R k > 0. 


(1.3) 


One often uses the electroneutrality conditions on the boundary concentrations because the 
solutions are made from electroneutral solid salts, 

n n 

^2 z s L s = z sR s = 0. (1.4) 

S = 1 S=1 

The electrochemical potential fi k (X) for the k- th ion species consists of the ideal component 
fi]f(X) and the excess component / u^(A"), i.e., fi k (X) = fi l £(X) + ^(X). The excess elec¬ 
trochemical potential //jp( X ) accounts for the hnite size effect of ions. It is needed whenever 
concentrations exceed say 50 mM, as they almost always do in technological and biological sit¬ 
uations. The classical PNP model only deals with the ideal component /^(X), which reflects 
the collision between ion particles and water molecules and ignores the size of ions; that is, 

li k {X) = /4 d (.X) = z k e 0 <S>(X) + k B Tln (1.5) 

Lo 

where Cq is a characteristic concentration of the problems, and one may consider, 


C 0 = max {L k ,R k , sup |Q(X)|}. 

l<k<n A'e[ 0 ,i] 


( 1 . 6 ) 


For given V, Q(X), L k s and R k s, if (4>(A'), C k (X), J k ) is a solution of the boundary value 


problem (BVP) of (1.2) and (1.3), then the electric current X is 

n 

X — Gq ^ ^ Z s Cfs- 


(1.7) 


S=1 


For an analysis of the boundary value problem (1.2) and (1.3), we work on a dimensionless 
form. Set 


V o = max { sup V k (X)} and e r = sup £>(X). 

xe[o ,L] 


l<k<n Xg [ 0 ^] 


Let 


^2 _ k B T 

' el^Co ’ 

Q(X) 


Q{x) = 


Co 


tr{x) = XXl, *-*, = aw-MQ, 

£r L L l LJ 0 

= ^-^(X), c k (x) = ^P~, jl k = -L-ijl k , J k=j ^. 
k B T Co k B T LCqVo 


( 1 . 8 ) 
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In terms of the new variables, the BVP of (1.2) and (1.3) become, for k = 1, 2, • • • , n, 

^ 2 d 


£ 


h(x) dx 


fi r (x)h(x)-^(pj = - ^ ~2z s c s - Q(x), 

' * s=l 

= 0 , ~ J k = h { x ) D k{x)ck^rfa, 
dx dx 


(1.9) 


with the boundary conditions 


</(0) = V=-^V, C fc (0 ) = l k = I ± ] 0(1) = 0, c fc (l ) = r k = ^. 

KB-L C/0 Oo 


( 1 . 10 ) 


Remark 1.1. The actual dimensional forms of quantities have been used for all figures throughout 
the paper, i.e., for k = 1,2, 


Ck =CoCk (M = Molar), 
k sT 2 . 3 


Q = C 0 Q (M), 
k B T , 


$ x 10 3 (mV), V = V (mV), 

eo eo 

Mfc = (eo$ + hBTln(C fc /C' 0 )) X 10 21 (zJ) = </(x) + lnc fc (in fc s T), 
J k =LC 0 V 0 J k x 10~ 6 (MTz), X = LCoV 0 e 0 I x 10 12 (pA), 


( 1 . 11 ) 


where we take Co = 10M, X = 2.5ran, and V q = 2.032 x 10 9 m 2 /s, and, for diffusion constants 

am ), fork = 1,2, 

X\. =1.334 x 10 9 m 2 /s for Na + , or 

£> fc =2.032 x 10 ~ 9 m 2 /s for Cl~, or (1.12) 

£> fc =0.792 x 10“ 9 m 2 /s /or Ca 2+ . 

1.2 Setup of the problem. 

We give the precise setup of the problem considered in this work. More precisely, we assume 
(AO) The ionic mixture consists of two ion species with valences z\ = —z 2 = 1; 

(Al) Dk{x) = Hfc /or k = 1, 2 is a constant and e{x) = 1; 


(742/ Electroneutrality boundary conditions (1-4) hold; 

(A3) The permanent charge Q is piecewise constant with one nonzero region; that is, for a 
partition 0 < a < b < 1 of [ 0,1], 


Q{x) = 


Qi = Q3 = o, ie(o,a)u(M), 

Q 2 , x E (a, 6), 


(1.13) 


where Q 2 is a constant. 


We assume that e > 0 in system (1.9) is small. The assumption is reasonable since, if 
L = 2.5 nm(= 2.5 x 10~ 9 m) and Co = 10 M, then e « 10 -3 (|15)). The assumption that e is 
small enables one to treat the BVP (1.9) and (1.10) of the dimensionless problem as a singularly 
perturbed problem. 
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A geometric singular perturbation framework for analyzing BVP of cPNP systems was devel¬ 
oped first in [144139] for ionic mixtures with two types of ion species. This general framework was 
extended to arbitrary number of types of ion species successfully only when two special struc¬ 
tures of the PNP system were revealed ([3D])- One special structure is a complete set of integrals 
for the limit fast (or inner) system that allows a detailed analysis of singular layer component of 
the full problem. The other is a state-dependent scaling of the independent variable that turns 
the nonlinear limit slow (or outer) system to a linear system with constant coefficients. The co¬ 
efficients do depend on unknown flux densities to be determined as a part of the whole problem, 
and this is the mathematical reason for the rich dynamics of the problem. As a consequence of 
the framework, the existence, multiplicity, and spatial profiles of the singular orbits are reduced 
to a system of nonlinear algebraic equations that involves all relevant quantities altogether - the 
physical evidence why this framework works since “Everything interacts with everything else” in 
determining ion channel functions. This geometric framework with its extensions to include ion 
size effects to some extents in EDI ESI EH has produced a number of results that are central to 
ion channel properties ([29, EH EH HD 031 El!); for example, it was shown in [3T] that, in order 
to optimize effects of the permanent charge, the channel should have a short and narrow neck 
within which the permanent charge is confined; and, it was shown in [55 that, large permanent 
charge is responsible for the declining phenomenon - decreasing flux with increasing transmem¬ 
brane electrochemical potential. We refer the readers to the aforementioned papers for more 
details on geometric singular perturbation framework for PNP as well as concrete applications 
to ion channel problems. 

In this paper, we will apply some results and follow the notations in [13] and [33] for analytical 
results where the quantities are all in their dimensionless forms. Besides, for simplicity, we use 
the letters l, r and Q o where l\ = I 2 = l, r\ = r 2 = r, Q = 2Qo■ Also, for definiteness, we choose 
a = 1/3, b = 2/3 and h(x ) = 1. 


Remark 1.2. We remind the readers that the quantities V,l,r,Ck,Q,4>, fiki Jki and I are di¬ 
mensionless quantities corresponding to the dimensional quantities V, L, R,Ck, Q,&, Jki'Rk 
andl, respectively, obtained from (1.8). We switch from dimensional form to the dimensionless 
form and vice versa several times throughout the paper. 


2 Zero current problems with general diffusion constants. 

In this section, we study when and how fluxes and membrane potential produce current reversal. 
Throughout this section, in order to express the effects of diffusion constants on zero-current 
flux and reversal potential, we study and compare the results for different cases of diffusion 
constants where V\ = V 2 and where T> 1 2 , to indicate and emphasize the differences. 

Diffusion is the phenomenon through which solute particles spread or mix as a result of their 
potential energy. It is a spontaneous process that acts to eliminate differences in concentration 
and eventually leads a given mixture to a state of uniform composition. Fick’s first law [18] 

O r \2 

describes diffusion of uncharged particles by where c is the concentration, T> is the 

diffusion constant and t is time. Frequently, the determination of diffusion constants involves 
measuring sets of simultaneous values of t, c, and x. These measured values are then applied 
to a solution of Fick’s law to get the diffusion constants. Many techniques are available for 
the determination of diffusion constants of ions (charge particles) in aqueous solutions ([3] f!9, 
E, 35, 133] . etc.). When diffusion constants are equal, classical electrochemistry tells that many 
electrical phenomena “disappear” altogether, e.g., the “liquid junction” is zero. If the diffusion 
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constants of potassium and chloride are equal, classical electrochemistry says that KC1 acts 
nearly as an uncharged species. Indeed, that is the basis for the saturated KC1 salt bridge used 
in a broad range of electrochemical experiments for many years. Therefore, the equal diffusion 
constants case is quite degenerate. 

Experimental measurements are exclusively performed under isothermal conditions to avoid 
deviation of T> values. Nevertheless, even diffusion constants of certain ionic species may differ 
from one method to another, even when all other parameters are held constant. Many things 
becomes much more complicated mathematically when the diffusion constants are not equal, 
however. This complexity is what makes many biological and technological devices valuable. 
Some kinds of selectivity depend on the non-equality of diffusion constants, as well. 


Applying GSP theory to cPNP equations (1.1) for two ion species with diffusion constants 
k = 1,2, the authors of 


T> k , k = 1,2, the authors of [43] obtained an algebraic matching system with eleven equations 
and eleven unknowns for zero current problems and singular orbits on [0,1]. In order to be 
able to analyze it, they reduced the matching system to the case where two ion valences satisfy 
z\ = — z 2 . It follows from [33] that the reduced matching system for zero current / = J± — J 2 = 0 
when z\ = — Z 2 = 1 is, 


where 


G 1 {A,Q 0 ,9) = V and G 2 (A, Q 0 ,9) = 0, 


0,(A, «„,«) =«(in + In i) - (1 + 0 ) In ± + In S “ Q ° 

G 2 (A , Qo, 9) =9Qo hi 5a + 9Q .° - N. 


S b -Q o' 


( 2 . 1 ) 


( 2 . 2 ) 


S b + 9Q 0 


In above, 9 = , and A is the geometric mean of concentrations at x = a, that is, 


A = v/ci(a)c 2 (a). 


(2.3) 


Moreover, 


S a = y Qq + A 2 , S), = y Qq + B 2 , N = A — l + S a — S b . (2.4) 

In what follows, the numerical observations help us have a better understanding of the zero- 
current problems and complete the analytical results obtained in 


2.1 Zero-current flux. 

We aim to clarify the relationships of ion fluxes with permanent charge and diffusion constants 
when current is zero. Recall that fluxes J\ and J 2 are equal for this case and let J denotes it. 


For any permanent charge Q = 2Qo, once a solution ( A , V) of (2.1) is obtained, it follows from 
matching equations (see Appendix in [43]) that J is given by 


J = - 


6D\D 2 (A — l) QD\D 2 (r — B) 


{D\ + D 2 ) 


{Dr + D 2 ) 


(2.5) 


Sign of zero-current flux J . It was observed in [16] that the Nernst-Planck equation in (1.9) 
(with dimensionless quantities) gives, for k = 1, 2, 


a r 1 

D k Jo h{x)c k {x) 


z k V + In -. 

r 


( 2 . 6 ) 
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Therefore, the sign of flux Jdepends only on the boundary conditions l, r and V. Note that 
(2.6) holds for any condition, not just zero-current condition. 

For zero-current problem, V = V rev depends on l, r , D i, D 2 , and Q as well in general. So 
the sign of zero-current flux J seems to depend on all quantities and to be difficult to figure 
out. It is not the case. A consequence of (2.5) together with the results in Theorem 3.4 in 
is that: 


The zero-current flux J has the same sign as that of l — r. 

The latter follows directly from a statement in Theorem 3.4 in |43| that, for zero-current, 
l — A has the same sign as that oil — r. That is consistent with observations in Figure [l] where 
T> 1 = 1.334 x 10 -9 m 2 / sioi Na + is fixed, and V 2 varies from the same value (i.e., p = 1), to 
V 2 = 2.032 x 10 -9 m 2 /s for Cl - , (i.e., p > 1), to a random large value (i.e., p 3> 1). 




Figure 1: The function J = J{Q) for various values of p = D 2 /D 1 : The left panel for L = 2 mM 
and R = 5 mM; the right panel for L = 5 mM and R = 2 mM. 


Dependence of zero-current flux J on Q and TVs. Concerning the dependence of the zero- 
current flux J on Q, we have the following. 

(i) IfV 1 = V 2 , then the zero-current flux J is an even function in Q, and it is monotonically 
increasing for Q > 0. 


In this case, 0 = 0 and hence, it follows from (2.1)—(2.3) that A can be determined from IV = 0. 


Since N is an even function in Q 0 , A is an even function in Q 0 , and so is the zero-current flux 


J from (2.5) 


(ii) IfV 1 V 2 , then the zero-current flux J is not an even function in Q and the mono¬ 

tonicity of the zero-current flux J in Q seems to be more complicated. 


In this case, it can be seen that G 2 in (2.2) is not an even function in Q, and hence, the zero- 
current flux J is not. We would like to point out that, it follows from [54], for fixed p = T^/Th, 
no matter how large, one always has J —> 0 as Q —> ± 00 , that is consistent with the observations 
in Figure [l] 

(iii) Another fascinating result is that for small values of Q, the magnitude of p = V 2 /V 1 
affects the monotonicity of the zero-current flux J. 

In this case, if one fixes V\, and let V 2 increases from small values to Vo —> 00 , (i.e., p —> 00 ), 


then it follows from (2.5) that there is a meaningful change in the monotonicity of flux, for small 


values of Q, that is not intuitive. Let us consider the case where L < R and Q < 0 is small. It 




























follows from (2.1) and (2.2) that when Q increases, for the geometric mean of concentrations A, 
one obtains 

(a) A increases if p ~ 1 (that is 9 ~ 0), and consequently the zero-current flux J decreases; 

(b) A decreases if p 1 (that is 9 1), and hence, the zero-current flux J increases. 

Thus, depending on the size of p, the zero-current flux J may increase or decrease for Q < 0 

small, which is also consistent with the observations in Figure [lj The analysis for the case with 
L > R is similar. 


2.2 Reversal potential V rt 


Experimentalists have long identified reversal potential as an essential characteristic of ion chan¬ 
nels [23.;, 127]. Reversal potential is the potential at which the current reverses direction, i.e., 
V = <£(0) — 4>(L) that produces zero current T. Using dimensionless form of quantities (see 
Remark |1.2| ), it follows from (2.1) and (2.2), (where there are two ion species with valences 
z\ = — Z 2 = 1) that for general permanent charge Q = 2Q$ ^ 0 with arbitrary diffusion con¬ 
stants ([33]), the reversal potential is, 


t/■ a (\ S a + 6Qo l ^ ( A(Qq 1 9) S a — Qo 

= #( 1,1 +ln r ] ~ (1 + e) mfhJ) + ln sUo7 


(2.7) 


Range of Reversal potential V rev . For fixed l , r, and for any given Q = 2Qo, it follows from 
Theorem 4.2 in [33] that there exists a unique reversal potential V re v such that V rev < | In £1. 
As Q —> ±oo, then V rev gets close to the boundary values, i.e. V rev —> ±ln 

Zero reversal potential. One particular case is when the reversal potential is zero. To examine 
under what conditions one obtains V rev = 0, it follows from Theorem 4.2. in [43] that, 

(i) if D\ = -Do, then V rev (Q ) = 0 for Q = 0, 

(ii) if D\ < D 2 , then there is a Q < 0, such that V rev (Q) = 0, 

(hi) if Di > D 2 , then there is a Q > 0, such that V re v(Q ) = 0. 

Considering the second case in above, the observations in Figure [3] show that as p = D 2 /D 1 
increases, magnitude of the corresponding Q becomes larger. In fact, as p —> 00 , then Q —> — 00 . 

Reversal potential V rev (Q) for Q = 0. For Q = 0, one has V rev (jd) = 9 ln - from Theorem 4.2 
in [33] where 9 = (D 2 — Di)/(D 2 + D\). Therefore, 

(i) if Di = D 2 , then U re ,;(0) = 0, 

(ii) if Di / D 2 , then V rev (0) has the same sign as that of 9(1 — r). 

Let us consider the case where D\ < D 2 for a moment. In that case, V rev {t)) has the same sign 
as that of l — r. This is reasonable, since for V = 0 we have | Ji| < | J 2 I (since all but Jk/Dk 
are independent of Dk in ( |2.6[ )), and to help | Ji| more than | J 2 I to get Ji = J 2 for zero current 
condition, one needs to increase V when l > r (and decrease V when / < r), and that is why 
V r ev( 0) > 0 for l > r (and V re v(0) < 0 for l < r ). This is consistent with observations in Figure 
[3] as well. The analysis for the other case, i.e., when D\ > D 2 is similar. 

Monotonicity of V rev with respect to Q. It follows from Theorem 4.4 in [43] that, for 9Q > 
0, dgVrev has the same sign as that oil — r. This analytical result does not make conclusion 
about the case for 9Q < 0, though. The observations in Figures [2] and [3] show that result holds 
for any 9 and Q. Thus we have 

Conjecture: V rev is increasing in Q for l > r and decreasing in Q for l < r. 
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We remark that, in Figure [ 2 J we take L = 20 mM, R = 50 mM, and V\ = 1.334 x 10 _9 m 2 /s 
and T >2 = 2.032 x 10 _9 m 2 /s which are diffusion constants of Na + and Cl - respectively (see the 
solid line), and T>\ = 1.334 x 10 _9 m 2 /s and Vo = 0.792 x 10 _9 m 2 /s, where V 2 is the diffusion 
constants of Ca 2+ (see the dashed line). 



Figure 2: V 


Vrev(Q) decreases when L < R, independent of values of diffusion constants. 




Figure 3: The function V = V re v{Q) : 
panel for L = hmM and R = 2 mM. 


The left panel for L = 2 mM and R = 5m M; the right 


Dependence of V r 


on p = D 2 /D 1 . In terms of the dimensionless form of quantities (See 

on p = Z? 2 /D\ for effects of D\ and D 2 ■ It 


Remark 1.2), let us examine the dependence of V r 
follows from Proposition 4.6 in [33| that 


The reversal potential V rev is increasing in p if l > r and is decreasing in p if l < r. 

This feature reveals a fantastic aspect that is not intuitive immediately. Recall the equation 


(2.6). Given the boundary values and diffusion constants, the values one obtains for all terms 


in (2.6) except Jk are independent of D^ (@Q]). The relation surely holds for the zero-current 
condition, i.e., Ji = J 2 with V = V re v Now , let fix D\ and increase D 2 (so p is increasing). 
Then |^ 2 1 increases since all but in (2.6) are independent of D 2 . Consequently, to meet 
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zero-current condition, we need to increase | Ji|. Intuitively increasing V rev seems to lead to an 
increase in | J\\. This intuition agrees with the results only for l > r. In the other case, i.e., 
I < r, it is the exact opposite, though. That is, for l < r, it says, as p increases, V rev decreases. 
This counterintuitive behavior could be explained by the fact that ci(.x) depends on V rev , and 
reducing V rev could still increase |Ji|. In fact, l < r will result in reducing V rev , but c\{x) 
changes in a way that consequently increases | J\ \. 

To illustrate the counterintuitive behavior, we provide a n ume rical result in Figure [4j We 

Now, suppose that T>\ = 


1.1 


choose Co, L and T>\ for Na + , as it is mentioned in Remark 

0.792 x 10~ 9 m 2 /s, and consider the boundary concentrations L = 20 mM,R = 50mM and 
Q = 1 M. In this case V re „ = — 16.7657mU and J = —10.6184 MHz. Now, if we increase 
T >2 to T >2 = 2.032 x 10~ 9 m 2 /s which is Cl - diffusion constant, then V rev = —19.5527 mV and 
J = — 11.3146MF/T. These values make sense now, based on the above discussion. Note that, 
we just pictured the middle part of the channel in Figure [4] since the sides are almost identical. 
One should notice that it is hard to realize, from Figure [4| how L < R will result in reducing 
Vrev The complicated behavior discussed above convinces us that sometimes numerical results 
cannot necessarily help us if we do not know the whole truth from analytical results. 



Figure 4: Graphs of C±(X) when T>\ is fixed, but we increase V> 2 - 


2.3 A comparison with Goldman-Hodgkin-Katz equation for V rev - 

In this section, we first recall the GHK equation ([22, 23]), which relates the reversal potential 
with the boundary concentrations and the permeabilities of the membrane to the ions. If the 
membrane is permeable to only one ion, then that ion’s Nernst potential is the reversal potential 
at which the electrical and chemical driving forces balance. The GHK equation is a general¬ 
ization of the Nernst equation in which the membrane is permeable to more than just one ion. 
The derivation of GHK equation assumes that the electric field across the lipid membrane is 
constant (or equivalently, the electric potential 4>(x) is linear in x in the PNP model). Under 
the assumption, the I-V (current-voltage) relation is given by 


/ = 


Vj2 z l D k 

k =1 


rk - he ZkV 

1 - e z k v 
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For the case where n = 2 and z\ = — z-z = 1, the GHK equation for the reversal potential is 


V r G ev HK ( P ) = ^ 


r + pi 
l + pr' 


( 2 . 8 ) 


The assumption that the electric potential <f>(x) is linear is not correct when applied to chan¬ 
nels in proteins. That is because proteins have specialized structure and spatial distributions of 
permanent charge (acid and base side chains) and polarization (polar and nonpolar side chains). 
Experimental manipulations of the structure of channel proteins show that these properties con¬ 
trol the biological function of the channel. The GHK equation does not contain variables to 
describe any of these properties and so cannot account for the biological functions they control. 
A linear cj>(x) is widely believed to make sense without channel structure presumably, in par¬ 
ticular, Q o = 0. However, this is not correct, either. It follows from (2.7) for Qq = 0 that the 
zeroth order in e approximation of the reversal potential in this case is, 


Weu(0, p) = ^———r In -• (2.9) 

p + 1 r 

Figure [ 5 ] compares V rev {0, p) in 
seen, from the left panel that when l and r are close (for example L = CqI = 20 mM, R = Cqt = 
50 mM), then the two curves have almost the same behavior. However, when we reduce L from 
20 mM to 1 mM, then the right panel shows a significant difference between the two graphs. 


(2.9) with from the GHK-equation in (2 




Figure 5: V rev (Q = 0, p) vs V^ K {p): The left panel for L = 20mM and R = 50 mM; the right 
panel for L = ImM and R = 50 mM. 


In Figure [6j we arrange a simple numerical result for the case where Q / 0 to compare the 
graphs of V re v(Qi p), obtained from (2.7), for various values of permanent charge Q. We consider 
L = 20 mM, R = 50 mM, and 0 < p < 5 for some values of Q, i.e., Q = 0 M, 1 M, 10 M. 
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Figure 6: V rev (Q,p) with various values of permanent charges. 


2.4 Profiles of relevant physical quantities. 


It follows from (2.1) and (2.2) that for any given Q , once a solution (A, V) are determined, all 


the other unknowns can be determined. We consider the dimensional form of quantities, and 
fix (Q, L, R,'Di,T> 2 ) to numerically investigate the behavior of Cf(X) and 4>(X) throughout the 
channel. Figures [7] and [8] graph the cases with small permanent charge Q = 0.0001M = O.lmM 
when L = 20mM,R = 50 mM, T>\ = 1.334 x 10 _9 m 2 /s, and V >2 = 2.032 x 10 _9 m 2 /s. In this 
case, we obtain J = —72.7387MHz and V re v = —4.4820mP. 



Figure 7: The functions C^f X) (left) and ^(W) (right) with Q = O.lmM. 

Furthermore, Figures [9] and FigurcJlO] show graphs of concentrations, electrical potential, 
and electrochemical potentials versus X, where L = 20 mM, R = 50mM, Q = 2M, and diffusion 
constants are the same as previous one. In this case, we obtain J = —11.3146 MHz and 
V rev = -19.5527mW. 
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Figure 8: The functions Hi(X) and n 2 (X) are increasing for L < R and Q = 0.1 mM. 
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Figure 9: The functions C\ (X) and Cz(X) (left) and the function $(X) (right) for Q = 1M. 
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Figure 10: The functions H\{X) and ^(X) are increasing for L < R and Q = 1 M. 
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3 The current-voltage and current-permanent charge behaviors. 


Ionic movements across membranes lead to the generation of electrical currents. The current 
carried by ions can be examined through current-voltage relation or I-V curve. In such a case, 
voltage refers to the voltage across a membrane potential, and current is the flow of charged 
ions through channels in the membrane. Another important data is current-permanent charge , 
i.e. I-Q relation. Dependence of current on membrane potentials and permanent charge is 
investigated in this section for arbitrary values of diffusion constants. 

To derive the I-V and I-Q relations, we rely on P3] where the authors showed that the 
set of nonlinear algebraic equations is equivalent to one nonlinear equation depending on A, 


defined in (2.3). All other variables can be obtained from the special variable A. It is crucial to 


realize that this is a specific result, not possible in many cases. One can only imagine that the 
resulting simplification produces controllable and robust behavior that proved useful as evolution 
designed and refined protein channels. The reduction allowed by this composite variable can 
be postulated to be a “design principle” of channel construction, in technological (engineering) 
language, or an evolutionary adaptation, in biological language. In particular, the current I can 
be explicitly expressed in terms of boundary conditions, permanent charge, diffusion constants, 
and transmembrane potential in the special case that allows the definition of A. In what follows, 
we derive flux and current equations - when diffusion constants are involved as well - in terms 
of boundary concentrations, membrane potential, and permanent charge. The I-V, I-Q, J-V, 


and J-Q relations are investigated afterward in section 3.2 


3.1 Reduced flux and current equations. 


It was shown in m that the singular orbits of BVP ( |1.9| ) can be reduced to the algebraic 
equation 

_ « 

(3.1) 


r? In I 6 - — — N = 0, 


Sa-r] 

where B = l — A + r, and S a , Sb and N were defined in 


, and, 


_ ^ Qo /V , i - Qo) \ , 

“ nrjr V + r(S„ -Qa)l 


N 


ln f ■ 


(3.2) 


r{S a - Qo) 

Once A is solved from (3.1), we can obtain the flux densities and current equations as follows, 

for k = 1, 2, 


<4 :=J k {V,l,r,D 1 ,D 2 ) =3D k (l - A) (l + (—l) fc ^ 


Qo' 

I := /(V, l, r, D u D 2 ) =Ji -J 2 = 3 (l - A) (d, -D 2 -^-{D 1 + D 2 


(3.3) 


For any given (l,r, D\, D 2 ,Qq,V) there exists a solution for the flux J and current I. The 
numerical results in the next section give us more information on “current-voltage” and “current- 
permanent charge” relations. 


3.2 Current-voltage and current-permanent charge relations. 

Dependence of current on diffusion constants. Now we reveal a remarkable feature of the 
theoretical results that is not intuitive. Suppose that (/, r, Q) is given (V is still free!). It follows 


from (2.4) for the definition of N that there exists an A so that N = 0. It consequently follows 
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which 


from (3.2) that, if V = In 


implies, 


B(S a - Qq) 
A(S b - Qq) 


, then rj = 0. Therefore, from (3.3), I = 3(1 — A)(D\ — D 2 ), 


For some fixed parameters (l, r, V, Q), the sign of I depends on the sign of D\ — D 2 . 

Current-voltage relations or I-V curves. Figure 


11 


is a numerical simulation from (3.1) 
and (3.2) of the TV curves for some values of Q with T> 1 = 1.334 x 10~ 9 m 2 /s and T >2 = 


2.032 x 10 9 m 2 /s. One may conjecture, incorrectly based on the numerical observations, that 


the value of current X, obtained from (3.3), is unique for any V and Q. However, this is incorrect, 


in general, as mentioned in previous section. The existence of a unique solution is important 
since the opening and closing properties of channels might be thought to arise from non-unique 
solutions mm)- 



Figure 11: The function X = X(V) for L = 20 mM and R = 50 mM. 


Current-permanent charge relations or I-Q curves. Note that the equation (2.6), in di¬ 
mensional form is, 


Jk f 0 V k A(X)C k (X) dX ~ Mfc(0) A “ 1 ’ 2 - 

The sign of J b: is determined by the boundary conditions, independent of the permanent charge 
Q. Nevertheless, as seen in Figure [T2j the magnitudes of s, and consequently, the sign and 
the size of the current X do depend on the permanent charge Q in general. 

The other interesting counterintuitive observation for I-Q relations is the different ways 
current X behaves with respect to permanent charges Q, for various values of V. It follows 
from numerical observations (see Figure |f~2j for example) that there exists some V*(T>\,T> 2 ) > 0, 
independent of L, R, such that for any V where V > V*, there exists a unique maximum for the 
current, X max , that corresponds to V. Likewise, there exists some V*(Pi,X ) 2 ) < 0, independent 
of L, R, such that for any V where V < V*, there exists a unique minimum for the current, X m j n , 
that corresponds to V. 
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Q in M 


Q in M 


Figure 12: The function X = I(Q ) withV i < T> 2 -' The left panel for L = 20 mM and R = 30 mM; 
the right panel for L = 30 mM and R = 20 mM. 

Besides, we claim based on numerical observations (not proven though) that there exists 
V(X>i,X> 2 ) = min{V*, |V*|}, such that 

(i) for any given V where |V| > V, the corresponding current X is non-monotonic in Q, but 

(ii) for any V where |V| < V, the corresponding current X is monotonic in Q. 

In particular, it can be seen, in section |3.4| that current is monotonic in Q for V = 0. In the 
end, we would like to mention that the diffusion constants affect the values V* and 14 above. 


3.3 Flux-voltage and flux-permanent charge relations. 

Flux-voltage relations via J-V curves. Now, we study flux and membrane potential rela¬ 
tions, i.e., J-V curves, for given values of permanent charges Q. In all cells, maintenance of cell 
volume is essential for survival. The membrane potential is likely to be a key regulator of this 
process in many cells. The membrane potential, feeds into the cell volume, controls mechanism 
by changing the driving force for ionic fluxes. 

The flux J depends on concentrations, diffusion constants, and the electrical potential across 
the membrane. Without permanent charges, i.e., when Q = 0, the flux of one ion species is 
independent of the other quantities, based on the classical PNP models for dilute mixtures (as 
is well-known). With the presence of a permanent charge, as expected, the classical PNP model 
shows the dependence of the flux of one ion species on the other ion species m- In particular, 
the effects of permanent charges on ionic flows could be intricate, depending on the interactions 
between boundary conditions and the channel geometry. 

Flux-permanent charge relations via J-Q curves. Numerical observations (see Figure |l3j) 
express an interesting feature that is not easy to prove theoretically. We consider the case where 
l < r. One has, 

(i) if V < — In then J\ < 0, and if V > — In then J\ > 0; 

(ii) ifV< In l -, then J 2 > 0, and ifV > In then J 2 < 0. 

A similar discussion holds for the other case when l > r. 

We discuss the case (i) for J\ on the left panel in Figure [13} Note that if one takes l = 0.002 
and r = 0.003 in dimensionless forms (i.e., L = CqI = 20mV, R = C^r = 30 mV in dimensional 
forms), then In - = —0.4055. Hence, V = 0<—ln^ = 0.4055 in this case, and J\ < 0 as 
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stated in (i). However, if one considers V = 0.5 (i.e., in dimensional form V = 11.7696mF), 
then V = 0.5 > — In and J\ > 0. 




Figure 13: The functions J\ = J\ (Q) (the left panel), and J 2 = Ji(Q) (the right panel), with 
L = 20 mM and R = 30 mM. 


3.4 Zero-voltage current and zero-voltage fluxes. 


The different permeability of the membrane determines the zero membrane potential (voltage) 
to different types of ions, as well as the concentrations of the ions, the permanent charge, and 
the shape of the channel. The current I{V), and the fluxes Jk(V), for k = 1,2, are called 
zero-potential current and zero-potential fluxes respectively when V = 0. For any given value of 
membrane potential V, the formulas for the current I(V), for small and large values of permanent 
charge Q, are provided in m and |5'4'| respectively. It follows from m that for small values of Q, 
applying V = 0, zero-potential current I s (0), and zero-potential fluxes J|(0),(in dimensionless 
forms as mentioned in Remark |1. 2 [) are, 


i s ( 0) =(l 
Jki 0 ) =(i 


-r)(D 1 -D 2 ) 


-r) Dl + (- if 


(l — r ) 2 

(21 + r)(l + 2r) In 
3 (l ~ r) 2 P k 
2(21 + r)(l + 2r) In l - 


jQ + 0(Q"), 

r 

-Q + o(Q 2 ), 


k = 1,2. 


(3.4) 


Furthermore, for large positive values of Q = 2Qo, with v = ^ (where v is small), it follows 
from [54] that zero-potential current I l ( 0) and zero-potential fluxes ,/((()) are, 


I\ o)=- 




1 + r \ 2 n ^ 

—p - (/ — r)u 

y/l + 


+ 


3(1 


Tjk( ‘ r <Ui + .ff) + l {l2 “ r2) “ I 0 "' ^ lnr)( ' + r? Y + ° ( " 2) ’ 


(Vl + V (Vl + y/f) 

4(0) =6D 2 Vf M-' / ~ r) - W + QPi f lr(Vi-Vf) 

Yi + Vr (Vi + v^? x (Vl + v^) 


(3.5) 


+ X (l 2 - r 2 ) - |(In l - In r)(l + r) 2 jzz + 0(v 2 ). 
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It can be verified from equations (3.4) and (3.5), that for small and large values of Q, the 
zero-potential current 1(0) is increasing in Q when l < r and it is decreasing in Q if l > r. We 
recall again that the dimensionless quantities can be converted to the dimensional forms from 
(1.10) and Remark 0 

Figure 14 admits the above conclusion. Besides, it suggests that the monotonicity of 1(0) 
holds for all values of permanent charge, not only for small or large values. We emphasize that 
the monotonicity of current I with respect to permanent charge Q is just true for zero membrane 
potential, i.e., V = 0. Indeed, one should recall from section 3.2 that when b / 0, then the 
current I is not monotonic in Q. 




Figure 14: The function X = X(Q) for V = 0: The left panel for L = 20mM and R = 30 mM; 
the right panel for L = 30 mM and R = 20 mM. 


4 Conclusion. 

In this paper, we first recall the analytical results in [33] for arbitrary diffusion constants. To 
investigate the reversal potential problems, i.e., when the current is zero, we do numerical 
investigations based on the analytical results in [33], where many cases are studied analytically. 
We derive several remarkable properties of biological significance, from the analysis of these 
governing equations that hardly seem intuitive. 

Biophysicists are also interested in the relation of current-voltage (I-V), and current-permanent 
charge (I-Q), as well as reversal potentials problems. To do that, we first recall the analytical 
results in [13], for arbitrary diffusion constants, to drive the flux densities and current equa¬ 
tions explicitly. One way to characterize channels is the current at zero potential, that is when 
V = 0, that has practical advantages. Since it is usually easier to measure a large current than 
a vanishing one, we analyzed this case, as well. Furthermore, we briefly study the special cases 
of small and large permanent charge for zero voltage case, based on the analytical results of 
m and [S3] respectively. To bridge between small and large values of permanent charges, we 
numerically study I-V and I-Q relations for this case as well. 
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